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(57) Abstract: A method and system of providing segmentation of volumetric three-dimensional image data set such as MRA im- 
ages. Initially, a three-dimensional MRA volume is input An initial surface S is then generated by thresholding the input A signed 
distance function v to S is then generated, where v=v(x,t) and S is the zero level set of v(,0). The process proceeds by iteiatively 
updating v according to Formula (I) , the updating terrmnates at convergence or as deterrnined by an operator. S 4 is then denned to 
be the zero level set of the current distance function v' and reinitialize v' to be a distance function to S* The volume is continually 
iteranvcly updated such that a final distance function v is obtained. The first output obtained from this volume is a segmentation of 
vessels in the MRA data, obtained by computing the zero level set of v. 
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when smooth closed cuives in 2D were proven to shrink to a point under mean curvature 
motion. Mean curvature flow of any hypersuiface was framed as a level set problem.. Foi 
application to 

image segmentation, a vector field was induced on the embedding space, so that the 
5 evolution could be controlled by an image gradient field or other image data. The same 
results of existence, uniqueness, and stability of viscosity solutions were obtained for the 
modified evolution equations for the case of planar curves, and experiments on real-world 
images demonstrated the effectiveness of the approach 

Cuives evolving in the plane became surfaces evolving in space, called "minimal 

1 0 surfaces' 5 - Although the theorem on planar cuives shrinking to a point could not be extended 
to the case of surfaces evolving in 3D, the existence, uniqueness, and stability results of the 
level set formalism held analogously to the 2D case. Thus, the method was feasible for 
evolving both cuives in 2D and surfaces in 3D, Beyond elegant mathematics, spectacular 
results on real-world data sets established the method as an important segmentation tool in 

15 both domains- One fundamental limitation to these schemes has been that they describe only 
the flow of 

hypersurfaces, i e ., surfaces of co-dimension 1 

The problem of curve-shortening flow for 3D curves has been studied, and the level 
set technique has been generalized to arbitrary manifolds in aibitiary dimension Analogous 
2 0 r esults were provided and extend the level set evolution equation to account for an additional 
vector field induced on the space 



SUMMARY OF THE INVENTION 

Accordingly, the invention presents the first implementation of geodesic active 
2 5 contours in 3D. Specifically, the system and method of the invention use these techniques 
for automatic segmentation of blood vessels in MRA images. The dimension of the 
manifold is 1, and its co-dimension is 2 

The invention utilizes the fact that the underlying structur es in the image are indeed 
3D cuives and evolves an initial curve into the cuives in the data (the vessels).. In particular, 
30 the segmentation techniques of the invention are based on the concept of mean curvatur e 
flow, or curve-shortening flow, from the field of differential geometry . The proposed MRA 
segmentation method uses a mathematical modeling technique that is well-suited to the 
complicated curve-like structure of blood vessels The segmentation task is defined as an 
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energy minimization over all 3D curves and uses a level set method to search for a solution- 
The approach is an extension of previous level set segmentation techniques to higher co- 
dimension.. 

The invention thus sets foith a method of providing segmentation of volumetric 
5 three-dimensional image data such as MRA images Initially, a three-dimensional MRA 
volume is provided An initial surface S is then generated by thresholding input A signed 
distance function v to S is then generated, where v=v(x,t) and S is the zero level set of v( ,0). 
The process proceeds by iteratively updating v according to 

v, = /l(Vv(x,0,V-v(x,/)) +— Vv(x,0 -H' x — r , the updating terminates at convergence or as 

g |V/| 

1 0 determined by an operator.. S' is then defined to be the zero level set of the curr ent distance 
function v' and reinitialize V to be a distance function to SV. The volume is continually 
iteratively updated such that a final distance function v is obtained. The first output 
obtained from this volume is a segmentation of vessels in the MRA data, obtained by 
computing the zero level set of v. 

15 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG- 1 is a maximum intensity projection of a phase-contiast MRA image of blood 
vessels in a brain; 

FIG.. 2 illustrates level sets of an embedding function u, for a closed curve in R 2 ; 
20 FIG. 3 illustrates a single sequence showing eight successive stages of a tubular 

object undergoing mean curvature flow; 

FIG 4 A shows a curve having a tangent to C at p, the normal plane, the image-based 
vector, and its projection onto the normal plane; FIG. 4B shows a curve using the e-level 
set method; 

2 5 FIGs . 5A-C illustrate an evolving helix under mean curvature flow with additional 

vector field: target curve, initial level set, level set after evolution with endpoints 
constrained, respectively; 

FIG. 6 is an operational block diagram of a MR segmentation system utilizing the 
invention; 

30 FIG. 7 is a flowchart of the segmentation algorithm in accordance with the 

invention; 

FIG 8 is a flowchart showing the details of the surface and volume initialization 
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poition of the algorithm in accordance with the invention; 

FIGs 9A-D illustrate a vertical bar evolving into a segmentation of a first dataset; 

FIG. 10 is a single 3D dataset, the fust image in each row is the maximum intensity 
projection of the raw data, and the second image is the segmentation result fiom three 
5 orthogonal viewpoints; 

FIG. 1 1 illustrates an image of a partial segmentation of the first dataset in FIG.. 10, 
the colorscale is continuous from darkest to lightest intensities, with darkest indicating a 
radius of curvature < 1mm and lightest indicating a radius of curvature >2mm. 



10 DETAILED DESCRIPTION OF THE INVENTION 

Intuitively, mean curvature flow refers to some curve evolving in time so that at each 
point, the velocity vector normal to the curve is equal to the mean curvature vector This 
concept is normally defined for arbitrary generic surfaces, but only curves are necessary for 
the invention, so the definition has been restricted for purposes of illustration. 
1 5 More formally, let C t , t > 0 be a family of curves in R 2 or R 3 , N the normal for a 

given orientation. That is, C is a curve, and t represents the time parameter or the index into 
the family of curves, not position The mean curvature flow equation is then given by the 
vector equation 

C,=kN (1) 

2 0 with given initial curve C(0)=C o , K the curvature of the curve, and C t the time derivative of 

the curve It wili be appreciated that only ID curves are considered here, as opposed to 
evolving surfaces, the mean curvature is just the usual curvature of the curve.. This motion 
is also called "curve-shortening flow 5 ' since it is the solution, obtained by Euler -Lagrange 
equations, to the problem of minimizing curve length: 

25 ^JC(p)\dp 

where p is the spatial parameter of the curve 

The basic idea of the level set method is given to evolve a planar curve C . Define 
a function u: R 2 — » R so that C is a level-set of u. The convention that C is, in particular, the 
zero level set of u, although this choice is not necessary for the method. The function u is 

3 0 now an implicit r epresentation of the curve C. The advantages of this r epr esentation ar e that 

it is intrinsic (independent of parameterization) and that it is topologically flexible since 
different topologies of C are represented by the constant topology of u. Let C 0 be the initial 
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curve. 

Evolving C according to 

t>/W (2) 
with initial condition C( ,0) = C 0 ( ) for any function (J, is equivalent to evolving u according 
5 to 

u,=p\Vu\ (3) 

with initial condition u( ,0) = Uo(-) and u 0 (C 0 ) = 0.. 

This result is independent of the choice of function u As is customary in the 
literature, u to be the signed distance function to the curve C as shown in FIG.. 2.. FIG. 2 
1 0 illustrates level sets of an embedding function u, for a closed cuive in R 2 . 

The level set evolution equations that follow enable one to evolve space curves, with 
evolution driven by both mean curvature and image information. In the following 
discussion, C is a curve in 3D . 

Let v:R 3 -»[0,oo) be an auxiliaiy function whose zero level set is identically C, that 
15 is smooth near C, and such that Vv is non-zero outside C. Foi a nonzero vector qeR", define 

P =I~&- 

as the projector onto the plane normal to q Further define l(Vv(x,t),V 2 v(x,t)) as the 
smaller nonzero eigenvalue of P Vv V 2 v P Vv . The level set evolution equation is then 

v t =X(7v(x f t) i V a v(x,t)). 

20 That is, this evolution is equivalent to evolving C according to C,=kN in the sense that C is 

the zero level set of v throughout the evolution 

FIG 3 demonstrates this evolution by illustrating evolving curves under mean 

curvature flow. FIG 3 illustrates a single sequence showing eight successive stages of a 

tubular object undergoing curve-shortening flow (mean curvature flow), where the curve is 
2 5 the centerline of the tubular shape The bumps are first smoothed out until the shape 

appr oximates a torus, then the torus shrinks to a point. 

The situation where there is an underlying vector field driving the evolution, in 

combination with the curvature term will now be described. Assume the desired evolution 

equation is of the form 
30 C t =KN-nd 

where TI is the projection operator onto the normal space of C (which is a vector space of 
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dimension 2) and d is a given vectoi field in R 3 . The evolution equation for the embedding 
space then becomes 

v=A(Vv 5 V 2 v)+Vvd. 
For the case of ID structures in 3D images, it is desirable to minimize 



\g(\vi(C(p))\)\cxppp 



where C(p):[0,l]->R 3 is the ID curve, l:[0,a]x [0,b]x[0,c]->[0,oo) is the image, and 
g:[0,oo)-»R + is a strictly decreasing function such that g(r)->0 as r-*x>„ For the current 
implementation, g(r)=exp(-r) is used because it works well in practice .. Another common 
choice is g(|VI|)=l/[l+|V I| 2 - By computing the Euler-Lagrange equations, it is found that 
10 the curve evolution equation is 

C,=KN--^n{H^-) y (4) 

g )v/| 

where H is the Hessian of the intensity function , The second term in the above equation is 
illustrated in FIG.. 4A, which shows a curve having a tangent to C at p, the noimal plane, the 
image-based vector, and its projection onto the normal plane. FIG 4B is a curve using the 
15 e-level set method. That is, 

e' V/ 
S M 

so the equation for the embedding space is 

v, = A(Vv(x 5 /),V 2 v(^0) + -Vv(x,r) (5) 

g l V/ l 

Thus, the mean curvature flow and level set methods have now been used to segment 
20 ID structures in 3D. FIGs 5A-C show how underlying image information can attract the 
evolving tube. FIGs 5A-C illustrate an evolving helix under mean curvature flow with 
additional vector field: target curve, initial level set, level set after evolution with endpoints 
constrained, respectively The underlying volumetric image data is shown, as a maximum 
intensity projection in FIG,. 5A. This volume was generated by drawing a cosine curve in 
2 5 the volume, then smoothing with a Gaussian filter.. FIG 5B shows the initial curve, a helix . 
The r esult of the evolution is shown in FIG. SC.. 

Segmentation refers to the task of separating foreground regions of an image from 
background r egions The system and method for segmentation of vessels from MRA using 
the described level set method in accordance with the invention will now be described FIG 
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6 is an opemtional block diagram of a MR segmentation system 600 utilizing the invention . 
The system includes a conventional MR scanner 602 lunning in conjunction with a MR 
computer that generates and stores raw MRA image slices . The MR data is then segmented 
by a computer 606. A flowchart of the segmentation algorithm in accordance with the 
5 invention is shown in FIG- 7.. 

The invention produces 3D suiface models of blood vessels based on magnetic 
resonance angiography (MRA) data. A scenario in which one could use the invention is 
illustrated in FIG. 6. Consider a patient who is suspected to have abnormalities in blood 
vessels in the brain, such as an aneuryism or other . Or perhaps one would like to view the 

1 0 blood vessels for study purposes . 

The patient's head is imaged in a magnetic resonance scanner 602. The image 
produced is a three -dimensional image. This means that it is actually a stack of many (often 
60) two-dimensional images, each of which pictures a "slice" or cross-section of the 
anatomical region imaged and stored by MR computer 604 , In this example, the head is this 

15 region. 

The invention utilizes a computer 606 to generate a 3D surface model of the blood 
vessels in the head, based on this 3D image This surface model could be displayed and 
manipulated on a standard computer . The surface model can be viewed by clinicians, 
radiologists, and other persons The surface model is often preferable to the raw 3D image 

20 in the areas of ease of interpretation, ease of further measurements, incorporation with other 
anatomical information, and other ar eas . 

Blood vessels appear in MRA images as bright curve-like patterns that may be noisy 
and have gaps. The data is a stack of slices where most areas are dark, but vessels tend to 
be bright. This stack is collapsed into a single image for viewing by performing a projection 

25 through the stack that assigns to each pixel in the projection the brightest voxel over all 
slices.. 

The task of MRA segmentation is complicated by the presence of imaging artifacts 
which appear visually similar to true vessel structures and also to partial voluming, the case 
of a small imaging area having a br ightness value that is a combination of the brightness 
30 values of vessels and of background because the vessel is only partially inside the area 
imaged. 

This specific segmentation problem is part of the high-level problem of developing 
computerized techniques for- the analysis of medical images. Automatic and semi-automatic 
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techniques can potentially assist clinicians and radiologists, saving them much of the time 

required to manually segment large datasets, or mor e generally facilitating measurements 
and interpretation of the images.. 

The 8-level set method is now described. Since the projection operator P q is defined 
5 only for non-zero vectors q, the method is undefined at Vv=0, which is the curve itself, and 
is numerically unstable near the curve. For this reason, v is regarded as a distance function 
to a "tube" of small radius e around the curve, instead of extracting the true ID curve. That 
is, 8-level set is evolved instead of evolving the true curve as shown in FIG. 4B. It will be 
appreciated that e does not denote a fixed value here.. The evolving shape is a "tubular" 

1 0 surface of some (unspecified and variable) nonzero width In addition to being more robust, 
this method better captures the geometry of blood vessels, which have nonzero diameter 

Banding will now be described., Instead of evolving the entire volume, only 
the portion of the volume within a narrow band of the zer o level set is evolved (the current 
surface), Noimally, the band is set to include voxels that are up to 6 voxels away from the 

1 5 surface The advantage of this technique is efficiency, and the disadvantage is that structures 
that ate outside the band may be missed if the potential function g does not have a large 
enough capture range to attract the segmentation to these structures.. The interpretation of 
banding is different from that in previous level set methods; they propagate image 
information from the zero level set to the rest of the bands, while the invention uses image 

20 information at each point 

A description of the use of curvatur e instead of eigenvalues is now provided . For 
computational efficiency and because of numerical instability of the gradient computations 
and thus the evolution equation near Vv=0, the level sets of the function v flow in the 
direction of the normal with velocity equal to the sum of theii smaller principal curvature 

25 and the dot product of Vv with the image-based vector field d.. Therefore, the smaller 
curvature is computed directly from v instead of as an eigenvalue of P Vv V 2 v P^ v3 and the 
invention uses >=| Vv| p 2 , where p 2 is the smaller curvature . 

With reference now to FIGs 7 and 8, the algorithm used in the method of the 
invention is described 

30 Initially, the 3D MRA volume is loaded into the computer. The vessels appear 

bright, background appears dark. An initial surface S is generated either by thresholding the 
inputted data set or by using a previously generated surface . A signed distance function to 
S is then generated. This distance function is a 3D volume, v. The volume v is a function 
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of spatial location x . Since it will be repeatedly updated over time, it can also be considered 
a function of time That is, v=v(x,t) and S is the zero level set of v( ,0) . The detailed 
flowchart of the surface and volume initialization is shown in FIG.. 8 

The method continues by iteratively updating v according to Equation 5. The 
5 algoiithm may also incorporate the image scaling teim previously described and/or an 
orientation term. The process teiminates at convergence or as determined by the user. 

S is redefined periodically or when needed, to be the zero level set of the current 
distance function v' and reinitialize v' to be the distance function to S'. The process 
continues by updating the volume according to the previous step. 
1 0 The loop above yields a final distance function v. The first output obtained from this 

volume is a segmentation of the blood vessels, obtained by computing the zero level set of 
v- Second, the center lines are obtained as the local minima of the distance function. Third, 
estimates of vessel diameter are obtained as a by-product of the computation of X in 
Equation 5. 

15 To control the tr ade-off between fitting the surface to the image data and enforcing 

the smoothness constraint on the surface, an image scaling term imscale is added to Equation 
5 to obtain 

v, =A(Vv(x,0,V 2 v(jc,r))-f™sc^*^-Vv(x,0 (6) 

g |V/| 

imscale being set by the user or pie-set to a default value.. 
20 With respect to gradient directionality, because vessels appeal brighter than the 

backgr ound, the image term is weighted by the cosine of the angle between the normal to 
the surface and the gradient in the image.. This cosine is given by the dot pr oduct of the 
respective gradients of v and I, so the update equation becomes 

- Z(Vv(x. 0, V 2 v(jc, 0) + imscale * ( Vv V/) * £ Vv<>, t) H ~ . (7) 

g |V/| 

25 For example, if the two vectors point in the same direction, then the brighter region 

is inside the surface and the darker region is outside. The angle between the vectors is 0, 
whose cosine is 1, so the image term is fully counted.. However, if they point in opposite 
directions, the negative weighting prevents the evolving vessel walls from being attracted 
to image gradients that point in the opposite direction. 

30 As customary in level set segmentation methods, the volume v is periodically 

reinitialized to be a distance function.. The zero level set S is extracted, then each point in 
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the volume is set to be its distance to S. For the implementation of the invention, this 
reinitialization is itself a level set method. To obtain the positive distances, the surface is 
propagated outward at constant speed of 1, and the distance at each point is determined to 
be the time at which the surface crossed that point. A second step propagates the surface 
5 inward to obtain the negative distances analogously. With respect to the initial 

surface, FIG. 8 shows a flowchart of a more detailed portion of the algorithm used to 
generate the initial surface . This initial surface (and thus the initial volume) is normally 
generated by thresholding the MRA dataset. However , the method does not r equire that the 
initial surface be near the target surface but may use any initial surface. FIGs. 9A-D 
1 0 illustrate a vertical bar evolving into the segmentation of the first dataset in FIG 1 1 . 

As shown in FIG 8, the datasets may be pre-processed to reduce noise and smooth 
For the results presented here, the raw datasets were convolved with an isotropic Gaussian 
of a=0.5. 

The segmentations ar e post-processed to remove any surface patches whose surface 
15 area is less than some threshold (a parameter of the method) to eliminate patches 
corresponding to noise in the original data. 

With respect to vessel radii estimation, the lar ger principal curvature can be useful 
in measuring the radii of the vessels for a particular application, since radius is the inverse 
of curvature. Ihis curvature can be easily computed when the smaller principal curvature 

2 0 is computed for the segmentation.. The option to color-code the segmentations can be added 

based on vessel radii, as estimated from the local larger principal curvature of the tubular 
surface 

The results of the segmentation on four datasets in accor dance with the invention are 
demonstrated in FIG 10., FIG. 10 is a single 3D dataset, the fust image in each row is the 
25 maximum intensity projection of the r aw data, and the second image is the segmentation 
result from three orthogonal viewpoints. 

The following parameter s were used in these experiments . All settings were chosen 
empirically . For the method, imscale varied across the datasets depending on the noise 
present. A cleaning threshold c indicated the minimum surface area of connected 

3 0 components of the surface to be retained in the post-processing cleaning step 

FIG. 10 shows a single 3D data set (in maximum intensity projection) and 
segmentation result from three orthogonal viewpoints The result shown is obtained by 
thresholding the raw data set. 
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Finally, the capability to color-code the vasculature surface is demonstrated based 
on local curvature. With reference to PIG 1 1, it will be appreciated that for a ribbon-like 
vessel, the flatter sides shows a large radius, and the sharply cuived edges show a small 
radius,. In this image of a partial segmentation of the first dataset in FIG. 10, the color scale 
5 is continuous from darkest to lightest intensities, with darkest indicating a radius of 
curvature < lrnm and lightest indicating a radius of curvature £2mm.. The curvatures output 
by the evolution have been smoothed by a 3x3x3 filter prior to coloring the surface. 

Although the present invention has been shown and described with respect to 
several preferred embodiments thereof, various changes, omissions and additions to the 
1 0 form and detail thereof, may be made therein, without departing from the spirit and scope 
of the invention 

What is claimed is: 
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CLAIMS 



1 LA method of providing segmentation of volumetric three-dimensional image data 

2 set comprising: 

3 (a) inputting a volume of said volumetric three-dimensional image data set; 

4 (b) generating an initial surface S by thresholding said input volume; 

5 (c) geneiating a signed distance function v to S, where v=v(x,t) and S is the zero 

6 level set of v( ,0); 

7 & V/ 

7 (d)iteiatively updating v according to v, =A(Vv(a 5 /),V v(x,0) +— Vv(r,0 #7—7, 

g |v/| 

8 said updating teiminates at convergence or as deteimined by an operator; 

9 (e) defining S' to be the zero level set of the cunent distance function v' and 

1 0 reinitialize v' to be a distance function to S'; and 

11 (f) continuing updating the volume according to step (d) such that a final distance 

12 function v is obtained, wherein 

1 3 the first output obtained from said volume is a segmentation of vessels in the three- 

1 4 dimensional image data set, obtained by computing the zero level set of v. 

1 2 . The method of claim 1 , wherein the function g in the equation of step (d) can be 

2 any function such that g:[0 5 oo)~->R + is strictly decreasing and g(r)->0 as r->oo. 

1 3 The method of claim 1 , wherein an image scaling teim imscale is added to the 

2 equation of step (d) so that v, = A(Vv(x ,t)y 2 v(x 9 t)) + imscale Vv(x,/) Hr— . 

g |V/| 

1 4. The method of claim 3, wherein an orientation teim costerm = -Vv-VI, so that 

2 V| =MVv(x,t)y 2 v{xJ))-imscale*(Vv V/)*^-Vv(x,0 Ht=- v 

g i v/ l 



1 5. The method of claim 4 further comprising varying the frequency of reinitializing 

2 v to be the signed distance function to its zero level set 

1 6 . The method of claim 5, wherein any distance function algorithm can be used to 
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2 reinitialize v during the reinitialization step. 

1 7 The method of claim 6, wherein a predefined shape with a predefined distance 

2 function is chosen instead of thresholding the volumetric three-dimensional image data set. 

1 8. The method of claim 7 further comprising extracting the centeilines of the 

2 segmented volumetric three-dimensional image data set, said extraction is performed by: 

3 extracting the set of points that are inside the segmented surface and that are also 

4 local minima of the distance function v; and 

5 defining the hierarchical shortest path. 



1 9 The method of claim 8, wherein said step of defining the hierar chical shortest path 

2 comprises: 

3 (a) defining the set of primary nodes of a graph to be the extracted local minima and 

4 also secondary nodes to be the other voxels inside the vessels of the volumetric three- 

5 dimensional image data set; and 

6 (b) connecting nodes that are adjacent in the three-dimensional volume with the 

7 Euclidean distance between the nodes; 

8 (c) finding, for each primary node, the shortest path to another primary node; 

9 (d) combining connected paths into a single primary node and again finding the 

1 0 shortest path to another primary node; and 

1 1 (e) repeating steps (a)-(d) until entire structure is connected or no more connections 

12 are possible within the given segmentation. 



1 1 0 The method of claim 7 further comprising inferring local diameter information 

2 about the vessels of the volumetric three-dimensional image data set, said information being 

3 visualized by coloring the segmented surface according to local diameter, said diameter 

4 information being obtained directly in the algorithm since it is the inverse of the larger 

5 principal curvature of the surface and this curvature is found during the computation of X 

6 in the update rule for v. 



1 11. The method of claim 7 further comprising removing from the segmented surface 

2 any surface patches whose area or volume is beneath a predetermined threshold 
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1 12. The method of claim 8 further comprising removing from the segmented suiface 

2 any suiface patches whose area or volume is beneath a predeteimined threshold prior to 

3 extraction of the centerlines . 

1 13 The method of claim 9 further comprising removing from the segmented suiface 

2 any suiface patches whose area or volume is beneath a predeteimined threshold.. 

1 14.. The method of claim 1, wherein said volumetric thiee-dimensional image data 

2 set comprises MRA data. 

1 1 5. A system of providing segmentation of volumetric three-dimensional image data 

2 set comprising: 

3 (a) means foi inputting a volume of said volumetric three-dimensional image data 

4 set; 

5 (b) means for generating an initial suiface S by thresholding said input volume; 

6 (c) means for generating a signed distance function v to S, where v=v(x,t) and S is 

7 the zero level set of v(,0); 

8 (d) means foi iteiatively updating v according to 

? z' W 

9 v, = A(Vv(*,/) 3 V^v(x,/)) + — Vv(;c,/) Ht— r , said updating terminates at convergence or 

g |V7| 

1 0 as determined by an operator; 

1 1 (e) means for defining S' to be the zero level set of the curr ent distance function V 

1 2 and reinitialize v* to be a distance function to S'; and 

13 (f) means for continuing updating the volume accoiding to (d) such that a final 

1 4 distance function v is obtained, wherein 

15 the first output obtained from said volume is a segmentation of vessels in the thi ee- 

1 6 dimensional image data set, obtained by computing the zero level set of v. 



RNsnnnin <wn 



On7«JA1A1 I > 



WO 00/79481 



1 / 11 



PCT/US00/17282 




BNSDOCI& <WO 0079481 A 1 I > 



WO 00/79481 

it 



2 / 11 



PCT/US00/17282 



U 



U = 



« = 0 



nwQnnnirv ^wn 



Wp 00/79481 



3 / 11 



PCT/US00/17282 




BNSDOfiiD; <WO 



WO 00/79481 PCT/US00/17282 

4/11 




Wp 00/79481 



5 / 11 



PCT/US00/17282 




BNSDOCID: <WO 



007948 1A1 I > 



WO 00/79481 



6 / 11 



PCT/US00/17282 




BNSDOrir* <WD 0fl7Odfl1Al I > 



it I 



,i 

WO 00/79481 



7 / 11 



PCT/US00/17282 



MRA Image 



Generate 
Distance 
Function 



OR 



Reinitialize 

Distance 

Function 



Evolve 


V 


Extract 
Zero 

Level Set 







Segmentation 



Fia.i 



BNSnOCID: -cWO 



Wp 00/79481 



8 / 11 



PCT/USOO/17282 



OR 



MRA image 



Gaussian 
smoothing 



Threshold 



Generate 

pre-defined 

surface 



OR 



Distance 
function 



pre 



FIG. 8 



WO 00/79481 9 /H PCI/USOO/17282 




FIG. 1C F1C. ?£> 



BNSDOCIO <WO 0079481A1 I > 



WO 00/79481 

ii 



10 / 11 



PCT/US00/17282 




007Q4A1A1 i > 



WQ 00/79481 



11 / 11 



PCT/US00/17282 




BNSDOCID:<WO 0079481 A 1 1 > 



INTERNATIONAL SEARCH REPORT 



A.. CLASSIFICATION OF SUBJECT MATTER 

IPC 7 G06T5/00 



According to International Patent Classification (IPC) or to both national classification and IPC 



Into, tonal Application No 

PCT/US 00/17282 



5. FIELDS SEARCHED 



Minimum documentation searched (classification system tallowed by classification symbols) 

IPC 7 G06T 



Documentation searched other than minimum documentation to the extent that such documents are included In the fields searched 



Electronic data base consulted during the international search (name ot data base and where practical . search terms used) 

WPI Data, INSPEC, PAJ, EPO-Internal 



C DOCUMENTS CONSIDERED TO BE RELEVANT 



Category • 


Citation of document with indication where appropriate, of the relevant passages 


Relevant to daim No 


P,A 


US 5 920 319 A (AHN DAVID K ET AL) 


1-15 




6 July 1999 (1999-07-06) 




abstract 






column 2, line 55 -column 3, line 7 






column 3, line 52 - line 59 






column 5, line 40 - line 61 






column 6, line 41 - line 52 




P,A 


US 6 058 218 A (CLINE HARVEY ELLIS) 


1-15 




2 Hay 2000 (2000-05-02) 






abstract 






column 3, line 38 - line 41 






column 3, line 64 -column 4, line 35 






-/- 





)M Further documents are listed in the continuation of box C 



Patent family members are listed in annex 



• Special categories of cited documents : 

A' document defining the general stale of the art which is not 
considered to be of particular relevance 

' E' earlier document but published on or after the international 
filing date 

T" document which may throw doubts on priority daim(s) or 
which is cited to establish the publication date of another 
citation or other special reason (as specified) 

"O* document referring to an oral disclosure use. exhibition or 
other means 

"P" document published prior to the international filing date but 
I a ter than the priority date d aimed 



T' later document published after the international filing date 
or priority date and not in conflict with the application but 
cited to understand the principle or theory underlying the 
invention 

X* document of particular relevance: the claimed invention 
cannot be considered novel or cannot be considered to 
involve an inventive step when the document Is taken alone 

'Y* document of particular relevance; the claimed invention 
cannot be considered to involve an inventive step when the 
document is combined with one or more other such docu- 
ments, such combination being obvious to a person skiUed 
in (heart 

"&" document member of the same patent family 



Date of the actual completion of the international search 

3 November 2000 


Date of mailing of the international search report 

13/11/2000 


Name and mailing address of the ISA 

European Patent Office P B 5818 Patentiaan 2 
NL - 2280 HV Rijswijk 
Tel. f>31-70) 340-2040, Tx 31 651 epo nl 
Fax: (+31-70) 340-3016 


Authorized officer 

Gonzalez Ordonez, 0 



Form PCT/ISA/210 (second stmt) ( July 1 992) 



. <wo 



00784mA1 I > 



page 1 of 2 



INTERNATIONAL SEARCH REPORT 



Intt lonal Application No 

PCT/US 00/17282 



(^Continuation) DOCUMENTS CONSIDERED TO BE RELEVANT 



Category 0 Citation of document with indication where appropriate of the relevant passages 



Relevant to daim No 



VERDOCK B ET AL: "ACCURATE SEGMENTATION 
OF BLOOD VESSELS FROM 3D MEDICAL IMAGES" 
PROCEEDINGS OF THE INTERNATIONAL 
CONFERENCE ON IMAGE PROCESSING 
(ICIP),US,NEW YORK, IEEE, 

16 September 1996 (1996-09-16), pages 
311-314, XP000704040 
ISBN: 0-7803-3259-8 
page 312, paragraph 3 

KIM D K ET AL: "BOUNDARY SEGMENTATION AND 
MODEL FITTING USING AFFINE ACTIVE 
SURFACEMODEL" 

IEEETENCON - DIGITAL SIGNAL PROCESSING 
APPLICATIONS, US, NEW YORK, NY: IEEE, 
26 November 1996 (1996-11-26), pages 
147-150, XP000781671 
ISBN: 0-7803-3680-1 
abstract 

paragraph '0001! - paragraph '0003! 



1-15 



1-15 



MRI 



SNELL J W ET AL: "MODEL-BASED 
SEGMENTATION OF THE BRAIN FROM 3-D 
USING ACTIVE SURFACES- 
PROCEEDINGS OF THE NORTHEAST 
BIOENGINEERING CONFERENCE, US, NEW YORK, 
IEEE, 

vol. CONF. 19, 18 March 1993 (1993-03-18), 

pages 164-165, XP000399775 

abstract 

page 164, paragraphs 
MATHEMATICAL, FORMULATION 



1-15 



Form PCT/ISA/210 (continuation of second sheet) ( July 1 992) 



BNSDOCJD: <WO 0079481 A 1 1 > 



page 2 of 2 



INTERNATIONAL SEARCH REPORT 

• information on patent family members 



Inte Jonai Application No 

PCT/US 00/17282 



Patent document 
cited in search report 



Publication 
date 



Patent family 
members) 



Publication 
date 



US 5920319 



06-07-1999 



US 6058218 



02-05-2000 



us 


5782762 A 


21-07-1998 


AU 


6180498 A 


09-09-1998 


CP 

tr 


uyoiyyj a 


rt Q 1 *5 1 AAA 

08-1Z-1999 


wu 




i/-oo-iyyo 


AU 


706280 B 


10-06-1999 


AU 


4138096 A 


23-05-1996 


CA 


2202401 A 


09-05-1996 


EP 


0794728 A 


17-09-1997 


OP 


10507954 T 


04-08-1998 


WO 


9613207 A 


09-05-1996 


US 


6083162 A 


04-07-2000 


DE 


19851597 A 


12-05-1999 


JP 


11242739 A 


07-09-1999 



Form PCT/1SA/210 (paiani family annex* < July 1SS2) 
BNSDOCID: <WO 0079481 A 1 1 > 



(12) INTERNATIONAL APPLICATION PUBLISHED UNDER THE PATENT COOPERATION TREATY (PCT) 



CORRECTED VERSION 



(19) World Intellectual Property Organization 

I nler national Bureau 

(43) International Publication Date 
28 December 2000 (28,. 12.2000) 




PCT 



(10) International Publication Number 

WO 00/79481 Al 



(51) International Patent Classification 7 : G06T 5/00 

(21) International Application Number: PC.T/USOO/t 7282 

(22) Inter national f iling Date: 23 June 2(XX) (23 06 2000) 

(25) Filing Language: English 

(26) Publication Language: English 



(30) Piiorit> Data: 

607140 6W 



23 June 1999(23 06.1999) US 



(71) Applicants: MASSACHUSETTS INSTITUTE OF 
TECHNOLOC\ [US/US]; 77 Massachusetts Avenue.. 
Cambridge MA 021.39 (US) INST1TUT NATIONAL 
DE RECHERCHE EN INFORMAT1QUE ET EN 
AUTOMAT I QUE [FR/FR]: Domuine de Voluceau. Roc- 
quencoun. Boile postale 105. F-78153 Lechesnay (FR) 



ECOLE NATIONAL E DES PONTS ET CHAUSSEES 
[ — /FR]: - (FR) THE BRIGHAM AND WOMEN'S 
HOSPITAL, INC.. [US/US]: 75 Francis Street.. Boston 
MA 021 15 (US). 

(72) Inventors: LORIGO. Liana: 117 Sciarappa Street. #2. 
Cambridge. MA 02141 (US) CRIMSON, W., Eric, 
8 Hawihonie Road. Lexington. MA 02420-1713 
(US) FATJGERAS, Olivier: 969, avenue de Pierreieu 
F-06560 Valbonne (FR) KERIVEN, Rcnaud: 32. avenue 
Gilles F-94340 Joinville (FR). WESTIN, Carl-Fredrik: 
- KIKIN1S, Ronald: 20 Chapel Street, Apt#804A. 
Brookline. MA 02446 (US) 

(74) Agents: CONNORS, Matthew, E et al : Samuels. 
Gauthier & Stevens LLR 225 Franklin Street. Suite 3300. 
Boston MA 021 10 (US). 



(81) Designated State (national): JP 



[Continued on next pagt j 



^ (54) Title: MRA SEGMENTATION USING ACTIVE CONTOUR MODELS 



MRA Image 



I 



Generate 
Distance 
Function 



Evolve 


V 


Extract 
Zero 
Level Set 







OR 



Reinitialize 

Distance 

Fttnctkw 



Segmentation 



00 
Tf 

ON 

© 
O 



v, - -«V><».r>.V'v(».<»-> ^ v, (».'>"^| 



(I) 



(57) Abstract: A method and system of providing segmentation of volumetric three-dimensional image data set such as MRA im- 
ages Initially, a three-dimensional MRA volume is input. An initial surface S is then generated by thresholding the input A signed 
distance function \ to S is then generated where v=v(x t) and S is the zero level set of v( .0) The process proceeds by iteratively 
updating v according to Formula (I) the updating terminates at convergence or as determined hy an operator S is then defined to 
be the zero level set ol the current distance function v and reinitialize v" to be a distance function to S\ The volume is continually 
iteratively updated such that a final distance function \ is obtained. The first output obtained horn this volume is a segmentation of 
vessels in the MRA data obtained by computing the zero level set of v. 



BNSDOCID: <WO 0079481A1 IA> 



WO 00/79481 A 1 IIIUII1IH1IIIIIMIIIIIIU 



(84) Designated Slates (regional): European palem (AT. BE 

ch. cy.de dk. es. fi fr gb gr. ie. it. lu mc 

NL PT. SE) 
Published: 

— with inter national search tepori 

(48) Dale of publication of this corrected version: 

6 June 2002 



(15) Information about Correction: 

see PCT Gazette No 23/2002 oi 6 June 2002. Section II 

For tw o-letter code<> and other abbreviations refer to the "Guid- 
ance \otes on C odes and tbbrevtaiiom" appearing at tin begin- 
ning of each regttlai iswe oj the I'CJ Gazette 



6NSDOCID: <WO 00704«1A1 IA> 



WO 00/7948 1 PCT/US00/1 7282 

1 

MRA SEGMENTATION USING ACTIVE CONTOUR MODELS 



PRIORITY INFORMATION 

This application claims priority from provisional application Ser . No . 60/140,609 
5 filed June 23, 1999 and incorporated herein by reference in its entirety . 

BACKGROUND OF THE INVENTION 

The invention relates to the field of volumetric three-dimensional image data 
segmentation, and in paiticular to MRA segmentation.. 

10 Automatic and semi-automatic magnetic resonance angiography (MRA) 

segmentation techniques can potentially save radiologists large amounts of time requited for 
manual segmentation and can facilitate further data analysis. It is desirable to develop 
computer vision techniques for the segmentation of medical images. Specifically, the 
segmentation of volumetric vasculature images is considered, such as the magnetic 

1 5 resonance angiography (MRA) image pictured in FIG 1. 

As shown in FIG. I, blood vessels appear in MRA images as bright curve-like 
patterns which may be noisy and have gaps. What is shown is a "maximum intensity 
projection". The data is a stack of slices where most areas are dark, but vessels tend to be 
bright. Ihis stack is collapsed into a single image for viewing by performing a projection 

20 through the stack that assigns to each pixel in the projection the brightest voxel over all 
slices This image shows projections along three orthogonal axes- 
Thresholding is one possible approach to this segmentation problem and works 
adequately on the larger vessels The problem arises in detecting the small vessels. 
Thresholding cannot be used for the small vessels for several reasons. The voxels may have 

2 5 an intensity that is a combination of the intensities of vessels and background if the vessel 
is only partially inside the voxel. This sampling artifact is called partial voluming.. Other 
imaging conditions can cause some background areas to be as bright as other* vessel areas, 
complicating threshold selection. Finally, the images are often noisy, and methods using 
local contextual information can be more robust. 

30 Mean curvature evolution schemes for segmentation, implemented with level set 

methods, have become an important approach in computer vision. This approach uses 
partial differential equations to control the evolution The fundamental concepts from 
mathematics from which mean curvature schemes derive were explored several years earlier 
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when smooth closed curves in 2D were proven to shrink to a point under mean curvature 
motion. Mean cuivatuie flow of any hypeisurface was framed as a level set problem. For 
application to 

image segmentation, a vector field was induced on the embedding space, so that the 
5 evolution could be controlled by an image gradient field or other image data The same 
results of existence, uniqueness, and stability of viscosity solutions were obtained for the 
modified evolution equations for the case of planar curves, and experiments on real-world 
images demonstrated the effectiveness of the approach. 

Curves evolving in the plane became surfaces evolving in space, called "minimal 

1 0 surfaces". Although the theorem on planar curves shrinking to a point could not be extended 
to the case of surfaces evolving in 3D, the existence, uniqueness, and stability results of the 
level set formalism held analogously to the 2D case. Thus, the method was feasible for 
evolving both curves in 2D and surfaces in 3D Beyond elegant mathematics, spectacular 
results on real-world data sets established the method as an important segmentation tool in 

1 5 both domains.. One fundamental limitation to these schemes has been that they describe only 
the flow of 

hypersuifaces, i .e., surfaces of co-dimension L 

The pr oblem of curve-shortening flow for 3D curves has been studied, and the level 
set technique has been generalized to arbitrary manifolds in arbitrary dimension. Analogous 
2 0 results were provided and extend the level set evolution equation to account for an additional 
vector field induced on the space. 

SUMMARY OF THE INVENTION 

Accordingly, the invention presents the first implementation of geodesic active 
2 5 contours in 3D.. Specifically, the system and method of the invention use these techniques 
for automatic segmentation of blood vessels in MRA images The dimension of the 
manifold is 1, and its co-dimension is 2.. 

The invention utilizes the fact that the underlying structures in the image are indeed 
3D curves and evolves an initial curve into the curves in the data (the vessels) . In particular, 
30 the segmentation techniques of the invention are based on the concept of mean curvature 
flow, or curve-shortening flow, from the field of differential geometry. The proposed MRA 
segmentation method uses a mathematical modeling technique that is well-suited to the 
complicated curve-like structure of blood vessels The segmentation task is defined as an 
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energy minimization over all 3D curves and uses a level set method to search for a solution 
The approach is an extension of previous level set segmentation techniques to higher co- 
dimension. 

The invention thus sets foith a method of providing segmentation of volumetric 
5 three-dimensional image data such as MRA images.. Initially, a three-dimensional MRA 
volume is provided.. An initial surface S is then generated by thresholding input. A signed 
distance function v to S is then generated, where v=v(x,t) and S is the zero level set of v( ,0). 
The process proceeds by iteratively updating v according to 

v, = >1(Vv(jc s /),V v(x,0) + — Vv(jc>/) H: — r , the updating terminates at convergence or as 

g |V/| 

1 0 determined by an operator S' is then defined to be the zero level set of the current distance 
function v' and reinitialize v r to be a distance function to S\ The volume is continually 
iteratively updated such that a final distance function v is obtained. The first output 
obtained from this volume is a segmentation of vessels in the MRA data, obtained by 
computing the zero level set of v 

15 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG.. 1 is a maximum intensity projection of a phase-contrast MRA image of blood 
vessels in a brain; 

FIG 2 illustrates level sets of an embedding function u, for a closed curve in R 2 ; 
20 FIG. 3 illustrates a single sequence showing eight successive stages of a tubular 

object undergoing mean curvature flow; 

FIG. 4A shows a curve having a tangent to C at p, the normal plane, the image-based 
vector, and its projection onto the normal plane; FIG. 4B shows a curve using the e-level 
set method; 

2 5 FIGs.. 5A-C illustrate an evolving helix under mean curvature flow with additional 

vector field: target curve, initial level set, level set after evolution with endpoints 
constrained, respectively; 

FIG. 6 is an operational block diagram of a MR segmentation system utilizing the 
invention; 

30 FIG.. 7 is a flowchart of the segmentation algorithm in accordance with the 

invention; 

FIG. 8 is a flowchart showing the details of the surface and volume initialization 
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portion of the algorithm in accordance with the invention; 

FIGs. 9A-D illustrate a veitical bar evolving into a segmentation of a first dataset; 

FIG 10 is a single 3D dataset, the first image in each row is the maximum intensity 
projection of the raw data, and the second image is the segmentation result from three 
5 oithogonal viewpoints; 

FIG.. 1 1 illustrates an image of a partial segmentation of the first dataset in FIG. 1 0, 
the colorscale is continuous from darkest to lightest intensities, with darkest indicating a 
radius of curvature < 1mm and lightest indicating a radius of curvatur e >2mm. 

10 DETAILED DESCRIPTION OF THE INVENTION 

Intuitively, mean curvature flow refers to some curve evolving in time so that at each 
point, the velocity vector normal to the curve is equal to the mean curvature vector .. This 
concept is normally defined for arbitrary generic surfaces, but only curves are necessary for 
the invention, so the definition has been restricted for purposes of illustration . 
15 More formally, let C t , t > 0 be a family of curves in R 2 or R 3 , N the normal for a 

given orientation.. That is, C is a curve, and t represents the time parameter or the index into 
the family of curves, not position.. The mean curvatur e flow equation is then given by the 
vector equation 

C,=kN (1) 

2 0 with given initial cuive C(0)=C o , k the curvature of the curve, and C, the time derivative of 

the curve- It will be appreciated that only ID curves are considered here, as opposed to 
evolving surfaces, the mean curvature is just the usual curvature of the curve This motion 
is also called "curve-shortening flow 15 since it is the solution, obtained by Euler-Lagrange 
equations, to the problem of minimizing curve length: 

25 ^fiCXppp 

where p is the spatial parameter of the curve. 

The basic idea of the level set method is given to evolve a planar curve C Define 
a function u: R 2 — ► R so that C is a level-set of u. The convention that C is, in particular , the 
zero level set of u, although this choice is not necessary for the method.. The function u is 

3 0 now an implicit representation of the curve C . The advantages of this repr esentation are that 

it is intrinsic (independent of parameterization) and that it is topologically flexible since 
different topologies of C are represented by the constant topology of u Let C 0 be the initial 
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curve.. 

Evolving C according to 

C,=/W (2) 
with initial condition C( ,0) = C 0 ( ) for any function JJ, is equivalent to evolving u according 
5 to 

w,=/?M (3) 

with initial condition u( ,0) = u^ ) and Uo(C 0 ) = 0 

This result is independent of the choice of function u. As is customary in the 
literature, u to be the signed distance function to the curve C as shown in FIG 2.. PIG. 2 
10 illustr ates level sets of an embedding function u, for a closed cuive in R 2 . 

The level set evolution equations that follow enable one to evolve space curves, with 
evolution driven by both mean curvature and image information In the following 
discussion, C is a curve in 3D. 

Let v:R 3 -^[0,oo) be an auxiliary function whose zero level set is identically C, that 
1 5 is smooth near C, and such that Vv is non-zero outside C For' a nonzero vector qeR", define 

P-I-&- 
' tf 

as the projector onto the plane normal to q Further define X(Vv(x,t),V 2 v(x,t)) as the 
smaller nonzero eigenvalue of P Vv V 2 v P Vv The level set evolution equation is then 

v=^Vv(x,t),Vyx,t)>. 

2 0 That is, this evolution is equivalent to evolving C according to C =kN in the sense that C is 

the zero level set of v throughout the evolution.. 

FIG. 3 demonstrates this evolution by illustrating evolving curves under mean 

curvature flow. FIG 3 illustrates a single sequence showing eight successive stages of a 

tubular object undergoing curve-shortening flow (mean curvature flow), where the curve is 
25 the centerline of the tubular shape The bumps are first smoothed out until the shape 

approximates a torus, then the torus shrinks to a point . 

The situation where there is an underlying vector field driving the evolution, in 

combination with the curvature term will now be described.. Assume the desired evolution 

equation is of the form 
30 C,=KN-nd 

where IT is the projection operator onto the normal space of C (which is a vector space of 



BNSDOCID: <WO 0079431 A 1 1A> 



li 



WO 00/79,481 PCI/USOO/17282 

6 

dimension 2) and d is a given vector field in R 3 . The evolution equation for the embedding 
space then becomes 

v-Pi(Vv,V 2 v)+Vvd 
F oi the case of 1 D structures in 3D images, it is desirable to minimize 

5 )g(\VI(C(p))\)\C(p)ty 

o 

where C(p):[0,l]->R 3 is the ID curve, l:[0,a]x [0,b]x[0 ? c]->[0,oo) is the image, and 
g:[0,oo)->R + is a strictly decreasing function such that g(r)-»0 as i— >oo For the current 
implementation, g(r)=exp(-r) is used because it works well in practice Another common 
choice is g(|VI|)=l/[l+|V If . By computing the Euler-Lagrange equations, it is found that 
1 0 the curve evolution equation is 

c,=*^-^n(//S, (4) 
g |v/| 

where H is the Hessian of the intensity function.. The second term in the above equation is 
illustrated in FIG . 4A, which shows a curve having a tangent to C at p, the normal plane, the 
image-based vector, and its projection onto the normal plane. FIG 4B is a curve using the 
1 5 e-level set method , That is, 

S W 

so the equation for the embedding space is 

v, =A(Vv(x 5 0,V 2 v(x,0) + ^Vv(x,0 H^- { (5) 

g |v/[ 

Thus, the mean curvatur e flow and level set methods have now been used to segment 
20 ID structures in 3D . FIGs 5A-C show how underlying image information can attract the 
evolving tube. FIGs. 5A-C illustrate an evolving helix under mean curvature flow with 
additional vector field: target curve, initial level set, level set after evolution with endpoints 
constrained, respectively The underlying volumetric image data is shown, as a maximum 
intensity projection in FIG 5A. This volume was generated by drawing a cosine curve in 
2 5 the volume, then smoothing with a Gaussian filter .. FIG 5B shows the initial curve, a helix 
The result of the evolution is shown in FIG 5C. 

Segmentation r efers to the task of separating for eground regions of an image fr om 
background regions The system and method for segmentation of vessels from MRA using 
the described level set method in accordance with the invention will now be described FIG 
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6 is an opeiational block diagram of a MR segmentation system 600 utilizing the invention. 
The system includes a conventional MR scanner 602 running in conjunction with a MR 
computer that generates and stores raw MRA image slices.. The MR data is then segmented 
by a computer 606.. A flowchart of the segmentation algorithm in accordance with the 
5 invention is shown in FIG 7.. 

The invention produces 3D surface models of blood vessels based on magnetic 
resonance angiography (MRA) data. A scenario in which one could use the invention is 
illustrated in FIG 6. Consider a patient who is suspected to have abnormalities in blood 
vessels in the brain, such as an aneuryism or other. Or perhaps one would like to view the 

10 blood vessels for study purposes.. 

The patient's head is imaged in a magnetic resonance scanner 602. The image 
produced is a three-dimensional image. This means that it is actually a stack of many (often 
60) two-dimensional images, each of which pictures a "slice" or cross-section of the 
anatomical region imaged and stored by MR computer 604. In this example, the head is this 

15 region 

The invention utilizes a computer 606 to generate a 3D surface model of the blood 
vessels in the head, based on this 3D image, This surface model could be displayed and 
manipulated on a standard computer The surface model can be viewed by clinicians, 
radiologists, and other persons. The surface model is often preferable to the raw 3D image 

20 in the areas of ease of interpretation, ease of further measurements, incorporation with other 
anatomical information, and other areas 

Blood vessels appear in MRA images as bright curve-like patterns that may be noisy 
and have gaps . The data is a stack of slices where most areas are dark, but vessels tend to 
be bright This stack is collapsed into a single image for viewing by performing a projection 

2 5 through the stack that assigns to each pixel in the projection the brightest voxel over all 
slices. 

The task of MRA segmentation is complicated by the presence of imaging artifacts 
which appear visually similar to true vessel structur es and also to partial voluming, the case 
of a small imaging area having a brightness value that is a combination of the brightness 
30 values of vessels and of background because the vessel is only partially inside the area 
imaged.. 

This specific segmentation problem is part of the high-level problem of developing 
computerized techniques for the analysis of medical images. Automatic and semi-automatic 
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techniques can potentially assist clinicians and radiologists, saving them much of the time 
required to manually segment large datasets, or more generally facilitating measurements 
and interpretation of the images 

The e-levei set method is now described . Since the projection operator P q is defined 
5 only for non-zero vectors q, the method is undefined at Vv=0, which is the curve itself, and 
is numerically unstable near the curve . For this reason, v is regarded as a distance function 
to a "tube" of small radius e around the cuive, instead of extracting the true ID curve.. That 
is, E-level set is evolved instead of evolving the true curve as shown in FIG 4B. It will be 
appreciated that e does not denote a fixed value here. The evolving shape is a "tubular" 

1 0 surface of some (unspecified and variable) nonzero width. In addition to being more robust, 
this method better captures the geometry of blood vessels, which have nonzero diameter. 

Banding will now be described Instead of evolving the entire volume, only 
the portion of the volume within a narrow band of the zero level set is evolved (the curr ent 
surface) . Normally, the band is set to include voxels that are up to 6 voxels away from the 

1 5 surface. The advantage of this technique is efficiency, and the disadvantage is that structures 
that are outside the band may be missed if the potential function g does not have a large 
enough captur e range to attract the segmentation to these structures . The interpretation of 
banding is different from that in previous level set methods; they propagate image 
information from the zero level set to the rest of the bands, while the invention uses image 

2 0 informati on at each point 

A description of the use of curvature instead of eigenvalues is now provided For 
computational efficiency and because of numerical instability of the gradient computations 
and thus the evolution equation near Vv=0, the level sets of the function v flow in the 
dir ection of the normal with velocity equal to the sum of their smaller principal curvatur e 

25 and the dot product of Vv with the image-based vector field d. Therefore, the smaller 
curvature is computed directly from v instead of as an eigenvalue of P^W P Vv , and the 
invention uses X=\ Vvl p 2 , where p 2 is the smaller curvature. 

With reference now to FIGs 7 and 8, the algorithm used in the method of the 
invention is described . 

30 Initially, the 3D MRA volume is loaded into the computer The vessels appear 

bright, background appears dark. An initial surface S is generated either by thresholding the 
inputted data set or by using a previously generated surface.. A signed distance function to 
S is then generated. This distance function is a 3D volume, v. The volume v is a function 
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of spatial location x. Since it will be repeatedly updated over time, it can also be considered 
a function of time. That is, v=v(x,t) and S is the zero level set of v( ,0). The detailed 
flowchart of the surface and volume initialization is shown in FIG.. 8. 

The method continues by iteiatively updating v according to Equation 5. The 
5 algorithm may also incoiporate the image scaling teim previously desciibed and/or an 
orientation teim . The process terminates at convergence or as determined by the user. 

S is redefined periodically or when needed, to be the zero level set of the current 
distance function v' and reinitialize v' to be the distance function to S'. The process 
continues by updating the volume according to the previous step 
10 The loop above yields a final distance ftmction v. The fust output obtained from this 

volume is a segmentation of the blood vessels, obtained by computing the zero level set of 
v Second, the centerlines are obtained as the local minima of the distance function Third, 
estimates of vessel diameter are obtained as a by-product of the computation of X in 
Equation 5 

15 To control the trade-off between fitting the surface to the image data and enforcing 

the smoothness constraint on the surface, an image scaling term imscale is added to Equation 
5 to obtain 

v, = A(Vv(x,0, V 2 v(x,0) + imscale * ^ Vv(x,r) H ^ (6) 

g |V/j 

imscale being set by the user or pre-set to a default value. 
2 0 With respect to giadient directionality, because vessels appear* brighter than the 

backgr ound, the image term is weighted by the cosine of the angle between the normal to 
the surface and the giadient in the image.. This cosine is given by the dot product of the 
respective gradients of v and I, so the update equation becomes 

=A(Vv(xj),V 2 v(x,0) + iw$cflfe*(Vv.V/)»^Vv(x,0-//S. (7) 

g |W| 

25 For example, if the two vectois point in the same direction, then the brighter region 

is inside the surface and the darker region is outside . The angle between the vectois is 0, 
whose cosine is 1, so the image term is folly counted.. However, if they point in opposite 
directions, the negative weighting prevents the evolving vessel walls from being attracted 
to image gradients that point in the opposite direction.. 

30 As customaiy in level set segmentation methods, the volume v is periodically 

reinitialized to be a distance function. The zero level set S is extracted, then each point in 
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the volume is set to be its distance to S. For the implementation of the invention, this 
reinitialization is itself a level set method To obtain the positive distances, the surface is 
piopagated outward at constant speed of 1, and the distance at each point is determined to 
be the time at which the surface crossed that point A second step propagates the surf ace 
5 inward to obtain the negative distances analogously- With respect to the initial 

surface, FIG 8 shows a flowchart of a more detailed portion of the algorithm used to 
generate the initial surface. This initial surface (and thus the initial volume) is normally 
generated by thr esholding the MRA dataset.. However, the method does not require that the 
initial surface be near the target surface but may use any initial surface. FIGs, 9A-D 

1 0 illustrate a vertical bar evolving into the segmentation of the first dataset in FIG. 11. 

As shown in FIG 8, the datasets may be pre-processed to reduce noise and smooth. 
For* the results presented here, the raw datasets were convolved with an isotropic Gaussian 
of a=0.5.. 

The segmentations are post-processed to remove any surface patches whose surface 
15 area is less than some threshold (a parameter of the method) to eliminate patches 
corresponding to noise in the original data. 

With respect to vessel radii estimation, the larger principal curvature can be useful 
in measuring the radii of the vessels for a particular application, since r adius is the inverse 
of curvature. This curvature can be easily computed when the smaller principal curvature 

2 0 is computed for the segmentation . The option to color-code the segmentations can be added 

based on vessel radii, as estimated from the local larger principal curvature of the tubular 
surface 

The results of the segmentation on four datasets in accordance with the invention ar e 
demonstrated in FIG. 10. FIG.. 10 is a single 3D dataset, the first image in each row is the 
2 5 maximum intensity projection of the raw data, and the second image is the segmentation 
result from three orthogonal viewpoints. 

The following parameters were used in these experiments.. All settings were chosen 
empirically. For the method, imscale varied across the datasets depending on the noise 
present. A cleaning threshold c indicated the minimum surface area of connected 
30 components of the surface to be retained in the post-processing cleaning step. 

FIG. 10 shows a single 3D data set (in maximum intensity projection) and 
segmentation result from three orthogonal viewpoints. The result shown is obtained by 
thr esholding the r aw data set 
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Finally, the capability to color-code the vasculature surface is demonstiated based 
on local curvature With reference to FIG. 1 1, it will be appreciated that for a ribbon- like 
vessel, the flatter sides shows a large radius, and the sharply curved edges show a small 
radius.. In this image of a partial segmentation of the fust dataset in FIG, 10, the colorscale 
5 is continuous from darkest to lightest intensities, with darkest indicating a radius of 
curvature < 1mm and lightest indicating a radius of curvatur e >2mm„ The curvatures output 
by the evolution have been smoothed by a 3x3x3 filter prior to coloring the surface. 

Although the present invention has been shown and described with respect to 
several preferred embodiments thereof, various changes, omissions and additions to the 
1 0 form and detail thereof, may be made therein, without departing from the spirit and scope 
of the invention 

What is claimed is: 
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CLAIMS 



1 1 . A method of providing segmentation of volumetric three-dimensional image data 

2 set comprising: 

3 (a) inputting a volume of said volumetric three-dimensional image data set; 

4 (b) generating an initial surf iace S by thresholding said input volume; 

5 (c) generating a signed distance function v to S, where v=v(x,t) and S is the zero 

6 level set of v( ,0); 

£ l VI 

7 (d) iteiatively updating v according to v, = A(Vv(x,0* V 2 v(x 5 /)) + — Vv(x,/) H- — r , 

8 i v/ l 

8 said updating terminates at convergence or as determined by an operator; 

9 (e) defining S' to be the zero level set of the current distance function V and 

1 0 reinitialize v' to be a distance function to S 1 ; and 

1 1 (f) continuing updating the volume according to step (d) such that a final distance 

1 2 function v is obtained, wherein 

1 3 the first output obtained from said volume is a segmentation of vessels in the three- 

1 4 dimensional image data set, obtained by computing the zero level set of v . 

1 2 . The method of claim 1, wherein the function g in the equation of step (d) can be 

2 any function such that g:[0,co)-*R + is strictly decreasing and g(r)->0 as r-xx> 

1 3. The method of claim 1, wherein an image scaling term imscale is added to the 

2 equation of step (d) so that v, = A(Vv{x y t\V-v{x,t)) + imscale*— Vv(*,/) Hr— 

g jV/| 

1 4 The method of claim 3, wherein an orientation teim costerm = -Vv VI, so that 

£' V/ 



2 v, = /l(Vv(x,/), V v(x,0) ~ imscale * (Vv •• V/) * — Vv(x,f) H- — : 

8 N 



1 5.. The method of claim 4 further comprising varying the fr equency of reinitializing 

2 v to be the signed distance function to its zero level set. 

1 6 . The method of claim 5, wherein any distance function algorithm can be used to 
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2 reinitialize v during the reinitialization step 



1 7 The method of claim 6> wherein a predefined shape with a predefined distance 

2 function is chosen instead of thresholding the volumetric three-dimensional image data set 



1 8. The method of claim 7 further comprising extracting the centeilines of the 

2 segmented volumetric three-dimensional image data set, said extraction is performed by: 

3 extracting the set of points that are inside the segmented surface and that are also 

4 local minima of the distance function v; and 

5 defining the hierarchical shoitest path 



1 9 The method of claim 8, wherein said step of defining the hierar chical shoitest path 

2 comprises: 

3 (a) defining the set of primary nodes of a graph to be the extracted local minima and 

4 also secondary nodes to be the other voxels inside the vessels of the volumetric three- 

5 dimensional image data set; and 

6 (b) connecting nodes that are adjacent in the three-dimensional volume with the 

7 Euclidean distance between the nodes; 

8 (c) finding, for each primary node, the shoitest path to another primary node; 

9 (d) combining connected paths into a single primary node and again finding the 

1 0 shortest path to another primary node; and 

1 1 (e) repeating steps (a)-(d) until entire structure is connected or no more connections 

1 2 are possible within the given segmentation. 



1 10 , The method of claim 7 further comprising inferring local diameter information 

2 about the vessels of the volumetric three-dimensional image data set, said information being 

3 visualized by coloring the segmented surface according to local diameter, said diameter 

4 information being obtained directly in the algorithm since it is the inverse of the larger 

5 principal curvature of the surface and this curvature is found during the computation of X 

6 in the update lule for v 



1 1 1 The method of claim 7 further comprising removing from the segmented surface 

2 any surface patches whose area or volume is beneath a predetermined threshold 
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1 12 The method of claim 8 further compiising removing from the segmented surface 

2 any surface patches whose area or volume is beneath a predetermined threshold prior to 

3 extraction of the center lines. 

1 1 3 The method of claim 9 further comprising removing fiom the segmented surface 

2 any surface patches whose area or volume is beneath a predetermined threshold 

1 14- The method of claim 1 , wherein said volumetric three-dimensional image data 

2 set comprises MRA data- 

1 15. A system of providing segmentation of volumetric three-dimensional image data 

2 set comprising: 

3 (a) means for inputting a volume of said volumetric thr ee-dimensional image data 

4 set; 

5 (b) means for generating an initial surface S by thresholding said input volume; 

6 (c) means for generating a signed distance function v to S, where v=v(x,t) and S is 

7 the zero level set of v( ,0); 

8 (d) means for 1 iteratively updating v according to 

9 v f = A(Vv(x,/) > V^v(x,/)) + — Vv(x,/) — r , said updating terminates at convergence or 

g |V/| 

1 0 as determined by an operator ; 

1 1 (e) means for defining S' to be the zero level set of the current distance function v' 

12 and reinitialize v' to be a distance function to S'; and 

13 (f) means for continuing updating the volume according to (d) such that a final 

1 4 distance function v is obtained, wherein 

1 5 the first output obtained fiom said volume is a segmentation of vessels in the three- 

16 dimensional image data set, obtained by computing the zero level set of v. 



WO 00/79481 ' PCT/US00/17282 



1/9 




FIG, J 



FIG. 3 

SUBSTITUTE SHEET (RULE 26) 



BNSDOCin:<WO fW7fi4ftlAl IA> 



WO 00/79481 



2/3 



PCT/US00/17282 




FIG. 2 



SUBSTITUTE SHEET (RULE 26) 



" 4 



WO 00/79481 ' PCT/US00/17282 



3/9 




FIG. 4 A 




FIG. 4B 



SUBSTITUTE SHEET (RULE 26) 



WO 00/79481 " 



PCT/USOO/17282 



4/9 







. i illil 1 ■ 

. **r "V* 












• • .......... - - 



FIG. 5 A 
FIG. 5B 
FIG. 5C 



SUBSTITUTE SHEET (RULE 26) 



WO 00/79481 



PCT/USOO/17282 



5/9 




SUBSTITUTE SHEET (RULE 26) 



BNSDOCID: <WO nf>704«1A1 IA> 



WO 00/79481 * 



PCT/US00/17282 



6/9 



MR A IMAGE 



GENERATE 
DISTANCE 
FUNCTION 



EVOLVE 


V 


EXTRACT 
/ERG 
LEVEL SET 







OR j<- 



REINITIALIZE 
DISTANCE k" 
FUNCTION 



SEGMENTATION 



FIG. 7 



MR A IMAGE 



GAUSSIAN 
SMOOTHING 



OR 



THRESHOLD 



Sthr 



GENERATE 
PRE-DEFINED 
SURFACE 



OR 



DISTANCE 
FUNCTION 



J pre 



fig. 8 



SUBSTTTUTE SHEET (RULE 26) 



RN.9nnr.irv <wn (v>tcmaiai 



■% * 



WO 00/79481 



PCT/USOO/17282 



7/9 



FIG. 9 A FIG.9B 




FIG.9C FIG.9D 



SUBSTITUTE SHEET (RULE 26) 



t 



WO 00/79481' PCT/US00/17282 



8/9 




FIG. 10 



SUBSTITUTE SHEET (RULE 26) 



WO 00/79481 " 



PCT/US0O/17282 



9/9 




nrv7a.dBiAi i A»» 



SUBSTITUTE SHEET (RULE 26) 



i 

INTERNATIONAL SEARCH REPORT 



A CLASSIFICATION OF SUBJECT MATTER 

IPC 7 G06T5/00 



According to International Patent Classification (IPC) or to bom national classification and IPC 



Inti Jonal Application No 

PCT/US 00/17282 



B. FIELDS SEARCHED 



Minimum documentation searched (classification system followed by classification symbols) 

IPC 7 606T 



Documentation searched other than minimum documentation to the extent that such documents are included in the fields searched 



Electronic data base consulted during the international search (name of data base and where practical search terms used) 

WPI Data, INSPE'C, PAJ, EPO-Internal 



C DOCUMENTS CONSIDERED TO BE RELEVANT 



Category * 


Citation ol document with indication where appropriate of the relevant passages 


Relevant to claim No 


P.A 


US 5 920 319 A (AHN DAVID K ET AL) 


1-15 




6 July 1999 (1999-07-06) 






abstract 






column 2, line 55 -column 3, line 7 






column 3, line 52 - line 59 






column 5, line 40 - line 61 






column 6, line 41 - line 52 




P,A 


US 6 058 218 A (CLINE HARVEY ELLIS) 


1-15 




2 May 2000 (2000-05-02) 






abstract 






column 3, line 38 - line 41 






column 3, line 64 -column 4, line 35 






-/-- 





m 



Further documents are listed in the continuation of box C 



Patent family members are listed in annex 



* Special categories of cited documents : 

"A" document defining the general state of the art which is not 
considered to be of particular relevance 

"E' earlier document but published on or after the international 
filing date 

f document which may throw doubts on priority daim(s) or 
which is cited to establish the publication date of another 
citation or other special reason (as specified) 

"O ' document referring to an oral disclosure use exhibition or 
other means 

*P* document published prior to the International filing date but 
later than the priority date claimed 



T* later document published after the international filing date 
or priority date and not in conflict with the application but 
cited to understand the principle or theory underlying the 
invention 

X" document of particular relevance; the claimed invention 
cannot be considered novel or cannot be considered to 
involve an inventive step when the document Is taken alone 

Y" document of particular relevance; the claimed invention 
cannot be considered to involve an inventive step when the 
document is combined with one or more other such docu- 
ments, such combination being obvious to a person skilled 
in the art 

"&* document member ol the same patent lamHy 



Date of the actual completion ol the international search 

3 November 2000 


Date of mailing of the international search report 

13/11/2000 


Name and mailing address of the ISA 

European Patent Office P B. 5818 Patentlaan 2 
NL - 2280 HV Rijswijk 
Tel. (+3 WO) 340-2040. Tx 31 651 epo nl, 
Fax: (+31-70) 340-3016 


Authorized officer 

Gonzalez Ordonez, 0 



Form PCT/1S A/210 (second *r»et) ( July 1992) 



page 1 of 2 



INTERNATIONAL SEARCH REPORT 



2 



C (Continuation) DOCUMENTS CONSIDERED TO BE RELEVANT 


Category • 


Citation of document with indication where appropriate, of the relevant passages 


Relevant I 


o claim No 


A 
n 


VFRHftflf R FT fll • "ArrilRATF ^FGMFMTATTPlM 

OF BLOOD VESSELS FROM 3D MEDICAL IMAGES" 
PROCEEDINGS OF THE INTERNATIONAL 
CONFERENCE ON IMAGE PROCESSING 
(ICIP),US,NEW YORK, IEEE, 
16 September 1996 (1996-09-16), pages 
311-314, XP00O70404O 
ISBN: 0-7803-3259-8 
page 312, paragraph 3 


1 _ 


.1 C 


A 


KIM D K ET AL: "BOUNDARY SEGMENTATION AND 

rlUULL rillllib UolNu Mrr lnJt MO live. 

SURFACEMODEL" 

lEEETENCON - DIGITAL SIGNAL PROCESSING 
APPLICATIONS, US, NEW YORK, NY: IEEE, 
26 November 1996 (1996-11-26), pages 
147-150, XP000781671 
ISBN: 0-7803-3680-1 
abstract 

paragraph '0001! - paragraph '0003! 


1- 


•15 


A 


SNELL J W ET AL: "MODEL-BASED 
SEGMENTATION OF THE BRAIN FROM 3-D MR I 
USING ACTIVE SURFACES" 
PROCEEDINGS OF THE NORTHEAST 
BIOENGINEERING CONFERENCE, US, NEW YORK, 
IEEE, 

vol . C0NF 19, 18 March 1993 (1993-03-18), 

pages 164-165, XP000399775 

abstract 

page 164, paragraphs 
MATHEMATICAL, FORMULATION 


1- 


•15 



Form PCT7ISA/210 (continuation ci sscord sheet) ( Jufy 1992) 



page 2 of 2 



Ink ional Application No 

PCT/US 00/17282 



INTERNATIONAL SEARCH REPORT 

information on patent family members 



Inte tonal Application No 

PCT/US 00/17282 



Patent document 
cited in search report 


Publication 
date 


Patent family 
member(s) 


Publication 
date 


lie* crtonom a 

US 5920319 A 


Oo-O 7-1999 


US 


5782762 


A 


A 1 rt7 1 rt A A 

21-07-1998 






AU 


6180498 


A 


Art Art 1 nrtrt 

09-09-1998 






EP 


0961993 


A 


AO 1 O 1 rt Art 

08-12-1999 






WO 


9837517 


A 


Z/-08-I99o 






AU 


706280 


B 


10-06-1999 






AU 


4138096 


A 


23-05-1996 






CA 


2202401 


A 


09-05-1996 






EP 


0794728 


A 


17-09-1997 






JP 


10507954 


T 


04-08-1998 






WO 


9613207 


A 


09-05-1996 






us 


6083162 


A 


04-07-2000 



US 6058218 A 02-05-2000 0E 19851597 A 12-05-1999 

JP 11242739 A 07-09-1999 



Form PCT/ISAtf 10 (patent family arme«) (July 1992) 



This Page is Inserted by IFW Indexing and Scanning 
Operations and is not part of the Official Record 

BEST AVAILABLE IMAGES 

Defective images within this document are accurate representations of the original 
documents submitted by the applicant. 

Defects in the images include but are not limited to the items checked: 

□ BLACK BORDERS 

□ IMAGE CUT OFF AT TOP, BOTTOM OR SIDES 
P'faded TEXT OR DRAWING 

□ BLURRED OR ILLEGIBLE TEXT OR DRAWING 

□ SKEWED/SLANTED IMAGES 

□ COLOR OR BLACK AND WHITE PHOTOGRAPHS 

□ GRAY SCALE DOCUMENTS 

A LINES OR MARKS ON ORIGINAL DOCUMENT 

□ REFERENCED) OR EXHIBIT(S) SUBMITTED ARE POOR QUALITY 

□ OTHER: 

IMAGES ARE BEST AVAILABLE COPY. 
As rescanning these documents will not correct the image 
problems checked, please do not report these problems to 
the IFW Image Problem Mailbox. 



